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Abstract 

We have developed a classical two- and three-body interaction potential to simulate the hydrox- 
ylated, natively oxidised Si surface in contact with water solutions, based on the combination and 
extension of the Stillinger- Weber potential and of a potential originally developed to simulate Si02 
polymorphs. The potential parameters are chosen to reproduce the structure, charge distribution, 
tensile surface stress and interactions with single water molecules of a natively oxidised Si surface 
model previously obtained by means of accurate density functional theory simulations. We have 
applied the potential to the case of hydrophilic silicon wafer bonding at room temperature, re- 
vealing maximum room temperature work of adhesion values for natively oxidised and amorphous 
silica surfaces of 97 mJ/m^ and 90 mJ/m^, respectively, at a water adsorption coverage of approx- 
imately 1 monolayer. The difference arises from the stronger interaction of the natively oxidised 
surface with liquid water, resulting in a higher heat of immersion (203 mJ/m^ vs. 166 mJ/m^), 
and may be explained in terms of the more pronounced water structuring close to the surface in 
alternating layers of larger and smaller density with respect to the liquid bulk. The computed 
force-displacement bonding curves may be a useful input for cohesive zone models where both the 
topographic details of the surfaces and the dependence of the attractive force on the initial surface 
separation and wetting can be taken into account. 



I. INTRODUCTION 

Direct wafer bonding has emerged as an important technology for sihcon-based microelec- 
tronic and micromechanical systems.-i^ In particular, the preparation of silicon-on-insulator 
devices takes advantage of the strong adhesion between oxidised and hydrated Si surfaces 
to bond together crystalline Si wafers without the need for adhesives or high pressures.-"^ 
The first step of the bonding process is the preparation of hydrophilic Si surfaces, in which 
a ~0.5 nm thin layer of native oxide covers the Si bulk and is terminated by chemisorbed 
hydroxyl groups and hydrogen-bonded water molecules. Afterwards, the surfaces are put in 
contact at room temperature by applying a small pressure in a localised region of the wafers. 
This triggers the propagation of a bonding front to the whole surface area,- in a process 
that is driven by the formation of a hydrogen-bonded water network trapped between the Si 
surfaces.- Finally, the bonded system is annealed at high temperature to induce condensa- 
tion of adjacent silanols on opposite surfaces and diffusion of trapped water molecules away 
from the interface, forming a Si02 layer between the Si wafers. 

Crucial for the formation of a defect-free bonded interface is complete bonding during 
the room temperature contacting step. To this end, the prepared surfaces need to be suf- 
ficiently fiat and clean and the non-covalent interactions between them sufficiently strong. 
The experimental evidence points towards a dependence of the interaction strength on the 
amount of water physisorbed on the surfaces,- which is approximately 1 monolayer (ML) at a 
relative humidity of about 70 %.- Moreover, under similar humidity conditions, the adhesion 
strength depends substantially on the details of the oxide structure covering the surfaces. 
The measured work of adhesion for natively oxidised surfaces with 1-2 ML of adsorbed water 
is approximately 100 mJ/m^.^ In the case of thicker oxide layers, obtained after thermal ox- 
idation of the Si surface, lower values in the range 60-85 mJ/m^ have been reported.™ These 
are consistent with the bonding energy associated with the closure of cracks in vitreous 
silica glass.- Here the measured work of adhesion is about 75 mJ/m^ for relative humidities 
above 20 %, but is found to rapidly decrease at lower humidity. Notably, equal adhesion 
strengths could be measured in experiments of crack opening or closing,^ pointing towards a 
reversible adhesion mechanism governed by hydrogen bonds and ruling out the formation of 
covalent bonds between the surfaces at room temperature. This is supported by a molecular 
dynamics simulation of amorphous silica wafers, which suggests a water-mediated adhesion 



mechanism in the absence of siloxane bridges across the bonding interface.— 

In this work, we perform classical molecular dynamics (MD) simulations to investigate 
the crucial role of trapped water during the room temperature stage of hydrophilic Si wafer 
bonding. In particular, we aim to investigate surfaces terminated by a native oxide layer, 
whose structure and composition have been previously determined by means of extensive first 
principles molecular dynamics (FPMD) simulations based on density functional theory.— li^ 
For this purpose, we develop a classical force field capable of reproducing the interactions 
between all chemical species present in the Si/SiOjj/water interface system previously ob- 
tained by FPMD. Our starting points are the well-known Stillinger- Weber potential,— which 
has been widely used to simulate Si bulk crystal structures, and a potential developed by 
Vashishta and co-workers to study the structure and energetics of Si02 polymorphs.— Here 
we combine these two potentials, extending them to take into account the full range of Si 
oxidation states present in the native oxide layer as well as the interactions between the 
oxide layer and liquid water. 

After describing the details of the potential development (Section [TTl) and the technical 
details of the MD simulations (Section IIIIj) . our investigation will start with the study of 
the structural and energetic details of the interfaces of both amorphous silica and natively 
oxidised Si with liquid water (Section IIVAO . Subsequently, results of our investigation 
into the hydrophilic bonding of pairs of silica wafers and natively oxidised Si wafers will 
be reported and compared (Section lIVBp . Finally, all of our main findings will be briefly 
discussed and summarised in Section |Vl 

II. FORCE FIELD DEVELOPMENT 

A. Reference model structure 

Under normal atmospheric conditions, the Si surface is passivated by a hydroxylated 
native oxide layer less than 1 nm in thickness. In Ref. [ill, the oxidation of the Si(OOl) 
surface in a dry environment has been investigated by means of FPMD simulations based 
on density functional theory (DFT). In good agreement with medium-energy ion scattering 
experiments,— it was found that oxygen spontaneously adsorbs onto the bare surface up to a 
coverage of 1.5 ML. At higher coverages, further oxide growth becomes limited by diffusion 



of O2 molecules to the reactive Si/SiO^. interface. In a subsequent work, water molecules 
were found to spontaneously adsorb and dissociate on the resulting oxide structure, creat- 
ing further reactive sites for reactions with oxygen molecules.— The oxidised and partially 
hydrated surface structure obtained in this series of FPMD simulations after the adsorption 
of seven oxygen and two water molecules is shown in Fig. [H The resulting concentration of 
4 hydroxyl groups on a 1.2 nm^ surface is consistent with values ranging between 2.6 and 
4.6 OH/nm^ measured experimentally on amorphous silica surfaces.— >i^ The 13 oxidised Si 
atoms in the final structure present variable formal valence charges ranging from Si'^ to 
Si^^, in ratios which are roughly consistent with the available experimental literature.— We 
have performed a MuUiken population analysis of this structure (projecting the ground state 
wave functions onto a basis consisting of s and p atomic orbitals) and found that the charges 
on the Si atoms increase approximately linearly with the number of nearest-neighbour O 
atoms, which act as electron acceptors (Table [T]). 

We have shown, in a recent work, that surface oxidation in a dry environment is accom- 
panied by the development of tensile surface stress up to 2.9 N/m at a coverage of 1.5 ML.— 
The surface stress g is defined as the basis invariant trace of the stress tensor a multiplied 
by half the height a of the periodically repeated supercell:— 

^ = -a Tr(a«/3) = -a{ayy + a^^), 



From the above equation, using the same DFT technique as in Ref. [181 (see Ref . I20I) . we have 
calculated the surface stress present in the relaxed hydroxylated native oxide, obtaining a 
value of 2.5 N/m (after subtracting the surface stress of the bottom hydrogen-terminated 
surface). This value is slightly lower than for the dry oxidised structure, since the dissociative 
adsorption of water on the native oxide has broken strained Si-0 bonds, although the surface 
clearly remains under tensile stress. 

B. Potential form for the natively oxidised Si surface 

We turn now to the development of a charge-based classical potential, whose analytic 
form and parameters will be constructed so as to reproduce the structure and tensile stress 
of the ah initio reference model described above. Our starting points are the well-known 
Stillinger- Weber potential^- (SW), which has been widely used to simulate bulk properties of 



silicon, and a potential developed by Vashishta et al. to simulate bulk and surface properties 
of silica polymorphsi^'-^^ (VP). Both of these potential energy functions contain only two- 
and three-body particle interactions and may be expressed by the general form: 

l<i<j l<i<j<k 

where the indices i and j run over all atoms of the simulated system and r^ is the vector 
connecting atom i with atom j. Our aim here is to combine the SW and VP potentials into 
a single form and extend them to include the full range of Si oxidation states present in the 
native oxide layer. In order to ensure that our potential reduces to the original forms^^*^ in 
the limits of bulk Si and bulk Si02, cut-off functions to the two- and three-body interactions 
will be introduced to smoothly interpolate between potential regimes. The form of these 
cut-off functions will be chosen so as to ensure that the potential energy and the forces 
on the particles are continuous functions of the particles' coordinates. This is essential to 
account for possible breaking and forming of Si-0 bonds in MD simulations of the surface 
with our potential form. 

In the VP potential, atomic charges with values of +1.6 e~ and —0.8 e~ are assigned to 
the Si^^ and O^^ species, respectively.— To maintain the original potential form in the case 
of Si^"*" species, we assign to the oxidised Si atoms and to the O atoms quantised partial 
charges according to: 



qsi = +0Aj2fcir,R„Ag) 
o 

qo = -0Aj2Mr,R„A,). 



Si 

According to the expression above, in agreement with the trend in computed Mulliken 
charges (see Table [T]), the charges on the oxidised Si species increase linearly with the number 
of oxygen neighbours up to a value of +1.6 for fully oxidised Si^+ species. The number of 
atomic neighbours is defined by the cut-off function fc, which falls off smoothly and with 
continuous derivative over a distance interval 2A around a cut-off distance R: 



fc{r,R,A) 



f 

1, r<R-A 

0, r>R + A. 



As anticipated above, the Si-Si two-body interactions are as in the SW potential, with 
the addition of a Coulomb term arising from the possible presence of charges on the Si 
atoms: 



Vsi-si' = MQsi, qsi')[Br ^ - l]exp(o-(r - a) ^)fc{r, Rsi, ^si) + 



QsiQsi 



r 

In the ab initio reference structure, the Si-Si bond lengths are observed to be approximately 
independent of the Si oxidation state. Therefore, in an adjustment to the original expres- 
sion,— the SW parameter A in the equation above includes a charge dependence, which 
counteracts the Coulomb repulsion between oxidised Si atoms (Fig. [2]): 

Aiq,q') = Ao{l + 3.2qq')fA{q)fA{q'). 

Si^^ species, as found in the tetrahedral form of Si02, form no direct Si-Si bonds and are 
instead connected via bridging oxygen atoms. Hence, the factor JaIq) is used here to remove 
the Si-Si attraction for Si^"*" species, recovering the purely repulsive behaviour of the original 
VP potential: 



fAiq) 



1.0 g<1.2 

4.0-2.5g 1.2 <q< 1.6 
0.0 g>1.6. 



The form of the Si-0 two-body interaction is based on Ref. 



14 and is made up of three 



terms: a short range steric repulsion, a charge-dipole attraction and a Coulomb attraction: 

Tr (CsioiQsi) Dsio I ,,n\ f j N , (lSi(lo 
ys^-o = (^ -^ ^exp(-r/6)J fBiqs^) + —^ ; 

Csio{.q) = Co-Ciq; 

, 2.5g q < 0.4 
fsiq) 

1.0 g>0.4. 

In the ab initio reference structure, the Si-0 bond lengths decrease approximately linearly 
with increasing oxidation state of Si. By varying the hard-core repulsion parameter Csio as 
a function of charge, we aim to reproduce the increase in bond length for oxidised Si atoms 
other than Si^"*" (Fig. [3]). The cut-off function /^(q') is introduced to reduce the uncharged, 
sub-surface Si atoms' interactions with O atoms in the oxide layer in order to recover the 
pure SW potential form within the Si bulk. 



The 0-0 interaction is exactly as in Ref. [ij and contains the same contributions as the 

Si-0 interaction: 

T/ ^oo' Doo' I /,x , <lo<lo' 
Vo-O' = — —eM-r/h) + —^ . 

The three-body interactions are all of the same SW form with parameters adapted from the 
literature,— >^>^ again smoothed by cut-off functions fc to recover the original VP potential: 

Vijk = Aexp[7i(rij - di)~^ + j2{rik - d2)~^]{cos9jik - cos9offc{rij, Rsi, ^si)fc{rik, Rsi, ^Si) ■ 

Finally, the hydroxyl groups that terminate the oxide surface are assigned charges equal to: 

qo = -0Aj2fcir,Rg, A,) -0.2 

Si 

Qh = +0.2 , 
consistent with the charge of approximately —0.6 assigned to the oxygen atoms of terminal 



silanol groups on the quartz surface in Ref. |23|. The 0-H bond length is constrained to 



0.975A and the Si-O-H interaction takes the above three-body form, with the parameters 
chosen to reproduce the variation about the minimum energy of the Si-O-H three-body 



interaction of Ref. 



23 



It should be emphasised that the potential form described here reduces exactly to the 
original SW and VP potentials for bulk Si and Si02 under equilibrium conditions. The 
cut-off distance Rsi is chosen to ensure that two- and three-body interactions involving Si 
are switched on when Si-Si separations are in the region of 2.6 to 2.8 A or below. In bulk 
Si (equilibrium bond length 2.36 A), our potential reproduces the SW potential up to bond 
lengths of 2.6 A, after which the attraction is smoothly reduced to zero, as shown by the 
dashed curve in Fig. [2l Similarly, our potential may be applied to Si02 polymorphs (Si- 
Si nearest neighbour peak separation of 3.1 A in amorphous silica, for example), reducing 
exactly to the VP potential for Si-Si separations greater than 2.8 A (apart from an additional 
Si-Si hard-core repulsion present in the original VP form, which is negligible for separations 
above 1 A). Below this Si-Si separation, we introduce repulsive Si-Si-0 and Si-Si-Si three- 
body interactions, while the Si-Si two-body interaction is removed entirely by the cut-off 
function /^(q'). 

Starting from the original values of the parameters for the separate SW and VP poten- 
tials,— >i^ the values of all parameters used in our two- and three-body interactions have been 
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carefully adjusted to reproduce as accurately as possible the ab initio reference structure in 
Fig. [1], and are listed in Table [TTl Average bond lengths and angles after minimisation of 
the hydroxylated surface structure, using this parameter set, show reasonable agreement 
with the ab initio minimised structure (Table IIIip . The calculated surface stress with our 
potential parameters is 2.0 N/m, consistent with the ab initio computed tensile stress of 
2.5 N/m. Importantly, no change in the topology of the SiO^: network (i.e. no breaking 
or forming of Si-0 bonds) was observed during 50 ps of classical MD simulation at 300 K, 
during which the charges and potential parameters were updated at every time step. 

As an additional test of our potential, we have performed a dynamical simulation of a 
Si/Si02 interface system composed of a 11 A thick slab of Si and a 17 A thick slab of a- 
quartz. At the interface, the Si(OOl) surface has been matched with the Si02 surface by 
expanding the a-quartz lattice parameter by 10 % in one direction to form an ideal, defect- 
free hetero junction. The system was separately annealed for 100 ps at 500 K, 100 ps at 700 K 
and for 50 ps at 1000 K. In none of these systems was breaking of existing bonds or formation 
of new bonds observed. Moreover, the average Si-0 bond lengths and Si-O-Si angles at 
the interface have been found to decrease from their bulk values by about 1 % and 4 % 
respectively, in agreement with a number of previous investigations.—"^ This demonstrates 
a certain degree of transferability of our potential to study Si/SiO^ heterogeneous systems. 
As a note of caution, we must say that we do not necessarily expect our potential to be 
predictive as far as the formation of unknown structures during chemical reactions (such as 
oxidation processes) is concerned. However, this may become possible upon augmentation 
of the potential within hybrid quantum/classical schemes such as the recently developed 
"Learn on the Fly" technique.— 

C. Interactions between the oxidised surface and water 

To study the behaviour of SiO^. surfaces in a wet environment, we model the interactions 
between the surface atoms and water molecules as a sum of Coulomb (via atomic partial 
charges) and Van der Waals (VdW) contributions. Consistent with force fields used to 
simulate solvated biomolecules in solution, the VdW interactions between water molecules 
and the surface hydroxyl groups are described by a "hydrogen bond" form V{rij) = A/rj^ — 
B/rj^, with the A and B parameters taken from standard AMBER parameter sets.^^ The 



VdW interactions between water molecules and all other atom types are of the standard 
Lennard- Jones (LJ) form: 

V{rij) = Aeij {{crij/nj^^ - (cTij/rijf) . 

For each pair of atoms we define: eij = Ju^j and Oij = 2^^/^(o"i+o"j), such that the minimum 
energy separation, tq, is (Xi + ctj, and V^ro) = —y/eiej. 

It has to be noted that the values of the partial charges on the surface atoms have been 
chosen so as to guarantee stability of the SiOa; network (see previous section). In particular, 
their values differ from the atomic point charges which best fit the electrostatic potential in 
a region outside the VdW radius of the atoms (ESP charges), which would be the best choice 
of charges to simulate the surface hydration properties.^ For this reason, we have adapted 
the parameters of the LJ part of the non-bonded surface-water potential so as to reproduce 
the binding energy curves for isolated water molecules on selected surface sites obtained in 
ab initio calculations. This energy is mostly accounted for by electrostatic effects, which are 
well described within standard DFT techniques. 

Within DFT,-^ we have calculated the total energy of a single water molecule as a function 
of its distance from the hydroxylated surface for two configurations, chosen such that the 
main interactions present are between one H atom of water with an O atom of the surface 
and between the O atom of water with a Si^"^ atom of the surface (Fig. H]). Identical 
total energy calculations were performed using the TIP3P water model^ and our newly 
developed force field for the native oxide, for different choices of the LJ parameters of the 
water-surface interactions. The best obtained binding energy curves are shown in Fig. HI 
and the corresponding, optimised set of LJ parameters is reported in Table IIVI As will 
be shown in Section IIVAI with this set of parameters the calculated heat of immersion of 
an amorphous silica surface amounts to 166 mJ/m^, which should be compared with the 
experimental value of 158 mJ/m^, measured for a surface with a density of terminal -OH 
groups of 3.4 nm~^.— 

III. COMPUTATIONAL DETAILS 

All classical simulations were performed using the DL_P0LY_3 MD package,— in which 
our newly-developed force field for the native oxide could be easily implemented. Long- 
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ranged Coulomb interactions were treated using the Ewald sum, with a real space cut-off of 
8 A. Newton's equations of motion were integrated with the velocity Verlet algorithm and 
a time step of 1 fs, using the RATTLE algorithn>2^ to constrain 0-H bond lengths. Where 
necessary, temperature control was achieved by velocity rescaling. 

An amorphous silica model was generated according to existing MD schemes,— using the 
VP potential in its original form,^^ which has been successfully applied to study amorphous 
silica surfaces.— Namely, we initially placed 1314 randomly positioned Si and O atoms, in 
the ratio 1 : 2, in a periodically repeated 30x30x20 A'^ simulation cell. The initial structure 
was annealed for 100 ps at 8000 K, then for a further 100 ps at 4000 K. Subsequently, the 
system was cooled to 300 K over a period of 360 ps, switching to a constant pressure ensemble 
controlled by a Nose-Hoover combined thermostat and barostat with 0.2 ps relaxation times. 
Following a further 200 ps annealing at 300 K, the bulk density was 2.37 g/cm^ (compared 
to the experimental value of 2.20 g/cm^).— Further structural analysis, detailed in Fig. [5] 
and Table |Vl was in good agreement with the literature^ and revealed the presence of a 
network of corner-sharing Si02 tetrahedra. 

Starting from the bulk amorphous silica structure obtained, a 50 A vacuum layer was 
inserted in the 2;-direction of the simulation cell, resulting in separated, periodically repeated 
slabs of surface area 9.1 nm^. The resulting surfaces were annealed at 1000 K for 500 ps, 
cooled to 300 K over 750 ps and, finally, annealed at 300 K for a further 500 ps, all at 
constant volume. These long annealing times were necessary to stabilise the concentrations 
of defects, found here in the form of threefold coordinated Si atoms and dangling O atoms, 
with converged concentrations of 1.6 nm~^ and 1.7 nm~^, respectively. These concentrations 



may be compared to Ref . |3J, which uses a Born-Mayer-Huggins pair interaction and a similar 
three-body potential, obtaining concentrations of the two defect types of 0.6 nm~^ and 
1.9 nm~^. Fig. [6l left, reveals that these defects were located in the first 5 A of the surface. 
In agreement with an existing VP potential simulation of the amorphous silica surface,— we 
found that the surface was terminated by dangling O atoms and that surface 0-Si-O bond 
angles were shifted to higher angles, indicating the presence of under-coordinated Si, as 
shown in Fig. [6|, right. After terminating these two defect types with -OH and -H groups, 
respectively, and relaxing the whole system, the resulting hydroxylated amorphous silica 
surface had a hydroxyl group surface concentration of 3.3 nm~^, which lies in the range of 
experimental measurements (2.6 - 4.6 nm~^).— li^ 
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To model the hydroxylated native oxide on Si(OOl), the ab initio reference structure shown 
in Fig. [1] was periodically repeated in the surface plane to form a 3 x 3 slab of surface area 
10.6 nm^. The bottom surface was terminated with a copy of the oxide layer after rotation 
through 180° and translation by a bulk Si lattice parameter. The resulting slab consisted of 
792 Si atoms, corresponding to eleven layers, oxidised and hydroxylated on either side with 
a total of 288 O and 72 H atoms. The lattice parameter was fixed to the equilibrium value 
for bulk Si obtained with the SW potential (5.44 A), and the vacuum gap in the direction 
perpendicular to the surface plane was set to 50 A, as in the case of the amorphous silica 
slab. 

To investigate water layering at the two hydroxylated surfaces (Section lIVAp . the 50 A 
vacuum layers were filled with TIP3P water molecules^ at a density of approximately 
1 g/cm^. For each surface, the entire system was equilibrated at 300 K and the height 
of the supercell adjusted to remove the stress in the direction perpendicular to the sur- 
face. Particle coordinates were collected every 100 time steps over a 1 ns production run. 
Selected snapshots from these simulations were used to construct the initial input files for 
the simulations of hydrophilic wafer bonding (Section II VB I) , keeping in the simulation cells 
only the water molecules closest to the surfaces up to heights corresponding to the desired 
surface water coverages. Based on the structural analysis of the surface/water interfaces 
(Section IIVAI) . one ML of adsorbed water is found to contain 98 and 117 water molecules 
on the amorphous silica and native oxide surface models, respectively, which correspond to 
surface densities of 10.7 and 11.0 molecules/nm^ . 



Based on Ref. 110|, force-displacement curves were calculated by reducing the height of 
the simulation cell perpendicular to the surfaces (and thus the separation between the top 
surface and the periodic image of the bottom surface) at a rate of 0.1 A every 11 ps. At each 
separation distance, the system was first equilibrated and thermalised to 300 K by velocity 
rescaling for 1 ps. During the subsequent 10 ps, the particle coordinates were collected every 
0.2 ps and the average z component of the stress in the simulation cell was calculated. In 
this way, the net interfacial force could be computed as a function of separation distance 
and the resulting curve integrated to obtain the total surface bonding energy. 
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IV. RESULTS AND DISCUSSION 

A. Structure and energetics of the SiO^; /Water interface 

We start our investigation of the interactions between wet oxidised silicon surfaces by 
studying the structure of bulk water in contact with both the amorphous silica and the 
natively oxidised Si surfaces. Any surface in contact with bulk water will have an effect on 
the intrinsic ordering of water in its proximity. In particular, water molecules will interact 
strongly with and partly penetrate into hydrophilic surfaces to form surface-water hydrogen 
bonds, in competition with water-water hydrogen bonds in the liquid. This will result in 
a layered structure of water molecules close to the surface, which can be quantified by the 
surface density profile of water molecules along the direction perpendicular to the surface. 

We have calculated the density profiles of water in contact with the two surface models by 
computing the average number of water molecules present in 0.1 A thin, planar slices parallel 
to the surface during a 1 ns simulation, as described in Section IIIII The results are reported 
in Fig. [3 In each case, we found evident structuring of water in layers of higher and lower 
density with respect to the bulk, consistent with previous simulations of a hydrophilic quartz 
surface.— For both surfaces, the position of the main peak in the density profile matches the 
position of the outermost surface hydroxyl group. The presence of subsidiary peaks closer to 
the surface, especially in the case of natively oxidised Si, is indicative of water penetration 
into the relatively open structure of the hydroxylated surfaces, as expected. With respect 
to the position of the first maximum (of density 1.05 g/cm^), the water density close to the 
amorphous silica surface presents a trough at a distance of 1.1 A and a second maximum at 
2.3 A, before reaching the bulk water density of 0.98 g/cm'^ at a distance of approximately 
4 A. The main peak at the native oxide surface is higher, reaching 1.32 g/cm^, and the 
oscillations continue out to approximately 5 A from the position of the main peak before 
reaching the bulk water density. 

The shape of the density profile obtained allows us to define a ML as the layer of water 
molecules contained between the surface and the position of the first trough after the main 
density peak. By integration of the density profiles, 98 and 117 water molecules are found 
to be contained in 1 ML on the amorphous silica and the natively oxidised surfaces, corre- 
sponding to surface densities of 10.7 and 11.0 water molecules/nm^ , or to surface areas per 
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water molecule of 9.3 and 9.1 A^, respectively. These values should be compared with early 
assumptions,- where an area of 10.6 A^ has been assigned to a water molecule adsorbed 
on quartz, the difference arising probably from the number of molecules which are able to 
penetrate relatively deeply into the surface oxide structures. 

For each surface, we have computed the heat of immersion as the difference between the 
energy of the system in contact with water and the energies of the two separate components, 
namely the dry surface model and bulk water. For each immersed system, the energy 
has been computed as the average value of the potential energy obtained in separate MD 
simulations at 300 K for 1 ns. With this approach, we calculated heat of immersion values 
of 166 mJ/m^ and 203 mJ/m^ for the amorphous silica and native oxide respectively. The 
first of the two values may be compared with experimentally measured values of 157 mJ/m^ 
or 158 mJ/m^, for surfaces with estimated densities of terminal -OH groups of 2.3 nm~^ or 
3.4 nm~^.— i^ The excellent agreement between the experimental and the theoretical values 
provides good support for the parametrisation of our surface-water potential, as described 
in Section III CI We note that a possible deprotonation of the surface hydroxyl groups in the 
experimental system is not expected to influence substantially the results obtained. In fact, 
in deionised water at pH 7, the surface charge is expected to be of the order of -1 /iC/cm^,— 
corresponding to just one deprotonated hydroxyl group per simulation cell, which is unlikely 
to contribute significantly to the computed heat of immersion. 

A further insight into the details of the surface/water interfaces may be gained by 
analysing the average number of surface-water and water-water hydrogen bonds formed 
per water molecule within the first ML compared with the liquid bulk. Considering a hydro- 
gen bond to be present between two O atoms when the 0-H- ■ ■ O angle is greater than 140° 
and the 0-0 separation is smaller than 3.5 A, the calculated average number of hydrogen 
bonds per molecule in bulk water is 3.13. Close to the surface, within the first ML, the 
corresponding values for the amorphous silica and the native oxide are 3.30 and 3.27, re- 
spectively, reflecting the competitive hydrogen bond formation at the surface. In each case, 
0.76 hydrogen bonds per molecule are formed with the surface (Table I VII) . which can be 
further divided between the bonds donated to O atoms bridging Si atoms of the surface, the 
bonds donated to O atoms of terminal hydroxyl groups, and the bonds received by H atoms 
of hydroxyl groups. The calculated values (Table IVII) indicate evident differences between 
the two surfaces. Namely, more hydrogen bonds are formed between water and bridging 
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oxygen atoms in the case of the native oxide structure, probably due to the reduced charge 
on some of the Si atoms, and thus the reduced electrostatic screening, compared with the 
amorphous silica surface. Since the charges on the bridging oxygen atoms are higher than 
those on the hydroxyl groups, this may help to explain the difference between the computed 
heat of immersion in the two cases, despite the small difference in the number of water 
molecules per surface area within the first ML. 

B. Wafer Bonding 

In the previous subsection, we have shown how competitive hydrogen bond formation 
at the surface/water interface leads to a local restructuring of the hydrogen bond network 
between the molecules. This results in oscillations in the water density, which propagate 
several Angstroms into the bulk liquid. The thin layer of water trapped between oxidised 
silicon wafers in the room temperature stage of hydrophilic wafer bonding will be subjected 
to the same effect, whose relationship to the attractive force between the wafers will be 
investigated below. 

We have calculated the average net force present between two bare amorphous silica 
surfaces and two bare native oxide surfaces in the absence of any water molecules as a 
function of the interface separation according to the procedure described in Section IIIII 
Starting from two separated surfaces, no attraction is observed when the surfaces are moved 
together (dashed lines in Fig. [8], top). When the surfaces become very close to each other, 
steric repulsion occurs between the surface atoms and the repulsive force increases rapidly. 
Zero separation in the curves in Fig. [8] was defined as the point of onset of this repulsive 
force. 

The situation is very different in the presence of water molecules between the surfaces. 
In Fig. [HI top, we have plotted the force-separation curves calculated at increasing water 
coverages (from approximately 0.25 to 2.5 ML per surface), along with the associated hy- 
drogen bond density as a function of surface separation (Fig. [HI bottom). As the surfaces 
move together, they experience a net attractive force (negative values in the reported plots), 
which initially increases, reaches a maximum value and then decreases, becoming repulsive 
at small separations. Subsequent debonding simulations, in which the surfaces were pulled 
apart at the same rate, revealed no hysteresis in the force-separation curves. 
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The overall behaviour of the system may be rationalised by looking at the number of 
hydrogen bonds formed per molecule between the surfaces. Starting from larger separations, 
the hydrogen bond density increases as the surfaces come into closer and closer contact. At 
smaller separations, if the concentration of water is smaller than about 1 ML (dotted lines 
in Fig. [S]), the interaction between the surfaces becomes repulsive well before the hydrogen 
bond density reaches the equilibrium values of 3.30 or 3.27 at the amorphous silica or native 
oxide surfaces, respectively (see Section HVAI) . In these cases, the computed force becomes 
dominated by the repulsion between the two solid surfaces, and the repulsive part of the 
force-separation curve tends to that calculated in the absence of trapped water. On the 
other hand, for water coverages greater than approximately 1 ML (solid lines in Fig. [S]), the 
onset of repulsion occurs before the surfaces interact directly. In fact, as is visible in Fig. [HI 
the point of zero force between the surfaces at increasing water coverages tends to the point 
at which the density of hydrogen bonds reaches 3.13, the equilibrium value in bulk water. 

These results imply that there is an optimum water concentration for room temperature 
hydrophilic wafer bonding, low enough that there is a high concentration of energetically 
favourable surface-water hydrogen bonds, as observed within the first ML of the systems 
in Section IIVA| but high enough that the oxide surface repulsion does not dominate. By 
integration of the force-displacement curves we can compute the energy gained during the 
simulated bonding process, that is the work of adhesion of the wet surfaces Qad- This is 
reported as a function of the water coverage in Fig. [9l indeed showing an optimum water 
concentration for bonding at approximately 1 ML for both surface models. At very low 
coverages, the steric repulsion between the surface dominates, so that Qad rapidly decreases 
to negligible values. At the other extreme, for high coverages, the values of Qad are very 
similar in the two cases and tend to the value of bulk TIP3P water, calculated by integration 
of a force-separation curve for a water-water interface, in the absence of any solid surface. 
The value obtained, 78 mJ/m^, is indicated with a dotted line in Fig. [9] and may be com- 
pared with the experimental surface energy of 72 mJ/m^ for the water-air interface at room 
temperature.— 

For the case of amorphous silica, the variation of Qad with water coverage may be com- 
pared with forces associated with crack closure in vitreous silica glass at different relative 
humidities, as reported in Ref. |6|. In this comparison, we shall take into account that Qad 
is half the strain energy release rate, and that each value of relative humidity is associated 
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with a well-defined coverage of adsorbed water.- The investigation in Ref. |6| shows that, at 
high humidity, Qad is roughly constant at about 75 mJ/m^, which corresponds to the surface 
energy of bulk water. At relative humidities lower than about 20 %, which corresponds to 
a water coverage of approximately 0.5 ML, Qad starts to decrease and drops to 25 mJ/m^ 
at a coverage of approximately 0.25 ML. This in good agreement with the shape of the 
amorphous silica curve shown in Fig. [9l where the computed work of adhesion is 19 mJ/m^ 
at 0.25 ML coverage, and is comparable with the surface energy of bulk water at large 
coverages. In our simulations, Qad presents a small peak of 90 mJ/m^ at approximately 
1 ML. Experimental values of work of adhesion in direct surface bonding experiments using 
silica glasses are lower, in the range of 60 to 85 mJ/m^,^ which might be explained by the 
microroughness present in the experimental wafer samples.— 

As far as the natively oxidised Si surface is concerned, our simulations indicate a slightly 
higher maximum value of Qad than amorphous silica, which may be related to the higher heat 
of immersion and to the more pronounced oscillations of the density close to the surface, 
as reported in Fig. [71 In fact, the water density profiles between closely spaced native 
oxide surfaces show oscillations that match those observed at an isolated slab (Fig. [TOl) . At 
0.5 ML coverage, the main peaks from opposite surfaces overlap, producing a single peak 
in the water density at equilibrium (Fig. [101 top right), fully consistent with trapping 1 
ML of water between the two surfaces. In this case, as shown in Fig. [TOl top left, the 
initial force required to separate the surface and thus break the hydrogen bond network 
of this ML is high, although Qad is relatively low, due to the surface- surface repulsion 
at small separation. At 1.0 ML coverage, Qad is maximum, reaching 97 mJ/m^, in good 
agreement with experimental values of approximately 100 mJ/m^ for natively oxidised Si 
wafers covered by 1-2 ML of water.-- In this case, the main peaks in the equilibrium water 
density at the two surfaces are separated by 2.5 A (Fig. [TOl bottom left), close to the ideal 
0-0 separation of 2.8 A for hydrogen-bonded TIP3P water.— This indicates the existence of 
an ordered, energetically stable, hydrogen-bonded water network, spanning the two bonded 
wafers. Finally, at higher coverages, the density in the central region between the wafer is 
constant, roughly at the equilibrium density of bulk water (Fig. [TOl bottom right), consistent 
with Qarf being comparable with the surface energy of liquid water. 

In an ideal debonding experiment, Qad may be thought of as being composed of a large 
energy loss (positive values) associated with breaking the network of hydrogen bonds between 
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the wafers and with an energy gain associated with the rearrangement of water at the 
surfaces, to optimise the number of hydrogen bonds per molecule. The more structured the 
water layer, the smaller will be the energy gain, due to the reduced capability of water to 
rearrange within the layer (which is consistent with the measured lower entropy of water 
molecules at low coverages on oxidised surfaces^). Therefore, the slightly larger Qad value 
than in the case of amorphous silica may be explained by the more pronounced structuring 
of water between the native oxide surfaces induced by the stronger surface-water interactions 
(see Section HVAI) . 



V. SUMMARY AND CONCLUSIONS 

In conclusion, we have developed a force field which enables us to investigate, at the 
classical level, the interactions between oxidised Si surfaces and a liquid environment. This 
is of particular importance given the increasingly wider fields of application of silicon-based 
microelectromechanical devices put in direct contact with water solution, such as sensors or 
actuators used in a physiological environment.-^ In our potential form, the surface interacts 
with the external environment through non-bonded interactions of the electrostatic and Van 
der Waals type, which makes the combination of the developed potential with standard force 
fields used to investigate biological macromolecules straightforward, and opens the way to 
large-scale simulations of biomolecule adsorption on realistic Si surfaces models. Moreover, 
the form of our force field allows processes associated with a rearrangement of the Si-0 
bond network to be simulated, as is necessary for the application of classical Hamiltonians in 
hybrid quantum-classical schemes of the "Learn on the Fly" type.— These appear to be very 
promising tools for the atomistic simulations of chemo-mechanical effects (for example, those 
occurring at a crack tip during a stress-corrosion process), which require both large system 
sizes and accurate quantum description of limited portions of the simulated systems.— "^"^ 

In this work, we have applied the newly developed potential to the simpler case of hy- 
drophilic Si wafer bonding at room temperature, where no breaking or forming of covalent 
bonds is expected to play a role. We have found that our parameterisation of the interatomic 
potential is able to reproduce the structural and mechanical properties of the natively oxi- 
dised Si surface, as obtained in previous FPMD simulations based on DFT.— i^ Moreover, 
calculated energetic details of the amorphous silica/ water interface, such as the surface heat 
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of immersion, were found to agree well with existing experimental measurements. As com- 
pared with the amorphous silica surface, the natively oxidised Si surface was found to interact 
more strongly with the liquid water, resulting in a higher heat of immersion (203 mJ/m^ vs. 
166 mJ/m^) and a more pronounced structuring of the water molecules close to the surface 
in alternating layers of larger and smaller density with respect to the liquid bulk. 

The layering of water at the oxidised surface has a profound effect on the computed 
force-displacement curves in simulated wafer bonding experiments at different coverages of 
adsorbed water. We have predicted that there is an optimum surface water concentration 
for bonding between both glassy silica and natively oxidised Si surfaces, for which the 
maximum work of adhesion values are 90 mJ/m^ and 97 mJ/m^, respectively, occurring 
at approximately 1 ML coverage (corresponding to about 70 % relative humidity-^). It must 
be noted that these values refer to perfectly flat surfaces, whereas experimentally measured 
values are expected to depend heavily on the nanoscale surface roughness present in the 
samples.— This is especially true in the case of bonding experiments, where a work of 
adhesion as low as 40 mJ/m^ has been measured for two natively oxidised Si surfaces.- The 
relatively slow diffusion of water parallel to the surface will prevent a perfect wetting of the 
whole surface during bonding. However, in subsequent fracture tests the work of adhesion 
has been shown to rise to 120 mJ/m^, possibly due to rearrangement of the water layer 
over the length scale of the surface roughness.— Our calculated force-displacement curves 
may provide a useful input for macroscopic models where both the topographic details of 
the bonded surfaces and the precise dependence of the attractive force on the initial surface 
separation and wetting can be taken into account. Combining the results presented here 
with these macroscopic models, such as the cohesive zone model of Kubair and Spearing,— 
will be the subject of future work.— 
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Tables 



TABLE I: Distribution of Mulliken charges in the reference structure of an oxidised and hydrated Si 
surface obtained in a series of FPMD simulations. Charges are multiples of the electronic charge. 



Species 


Mulliken Charges 


Average 


Si+ 


0.62 0.65 0.46 0.47 0.61 


0.56 


Si2+ 


1.26 1.28 1.28 


1.27 


Si3+ 


1.79 1.82 1.79 


1.80 


Si4+ 


2.30 2.34 


2.32 



TABLE II: Parameters of the interatomic potential. Units of length are A and those of energy are 
eVA(= 14.39 eV). 

Two-body potential parameters 



Rq 


2.0 


Ao 


1.062 


Co 


14.871 


Coo' 


51.692 


\ 


0.1 


B 


11.603 


Ci 


2.178 


Doo' 


1.536 


Rst 


2.7 


a 


3.771 


b 


4.43 






As^ 


0.1 


a 


2.0951 


Ds.o 


3.072 







Three-body potential parameters 



XsiSiSi 3.164 XsiSiO 3.164 Xsiosi 1-400 XosiO 0.350 XsiOH 0.521 

7i,SiSiSi 2.51 7i,5JSiO 4.06 ^i,siOSi 1-00 7i,osiO 1-00 ^i,siOH 1-00 

72,SiSiSi 2.51 ^2,SiSiO 0.52 ^2,SiOSi 1-00 j2,OSiO 1-00 ^2,SiOH 1-00 

di,SiSiSi 3.771 di^siSiO 3.981 di^siOSi 2.60 di^osiO 2.60 di^siOH 2.60 

d2,SiSiSi 3.771 d2,SiSiO 2.933 d2,SiOSi 2.60 d2,OSiO 2.60 d2,SiOH 2.60 

Oo,SiS^Si 109.47° eo,sis^o 109.47° eo,sios^ 141.00° ^0,05*0 109.47° Oo^siOH 122.50° 
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TABLE III: Comparison between ab initio and classical minimisations of the hydroxylated Si(OOl) 
surface. In the ab initio structure, Si-0 bond lengths (A) increase with decreasing Si oxida- 
tion state, while Si-Si bond lengths remain approximately constant. In the classical simulation, 
although some bond lengths and angles have changed, there is no restructuring of the surface. 
Si-O-H angles (126°) also agree with DFT values (127°). 





ab initio structure 




classical structure 




Si-0 Si-Si 0-Si-O Si-Si-0 Si-Si-Si 


Si-0 Si-Si 0-Si-O Si-Si-0 Si-Si-Si 


Si4+ 


1.64 - 109° 


- 


- 


1.63 - 109° 


- 


Si3+ 


1.65 2.36 109° 


109° 


- 


1.61 2.41 119° 


96° 


Si2+ 


1.67 2.39 109° 


108° 


109° 


1.62 2.44 139° 


101° 110° 


Si+ 


1.70 2.38 


108° 


111° 


1.65 2.39 


108° 111° 



TABLE IV: LJ parameters optimised for the surface-water interactions. Water OW parameters 
are taken from the TIP3P water model. Parameters for water HW atoms are for interactions with 
Si and OB only, otherwise ej is zero. In the surface, OB refers to bridging O atoms and OH to 
hydroxyl O atoms. 





water 


surface 




OW HW 


Si OB OH H 


ei/eVxlO-2 


0.66 0.13 


1.30 1.13 0.65 0.09 


CT./A 


1.768 0.650 


1.600 1.762 1.650 1.000 
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TABLE V: MD results for amorphous silica at 300K. The peak bond length or angle and the 
corresponding full width at half maximum (FWHM) of each distribution is reported and compared 



to Ref. 



14. 





this work 
bond length/A FWHM/A 


Ref. 14 
bond length/A FWHM/A 


Si-0 
0-0 

Si-Si 


1.62 0.07 

2.63 0.20 
3.09 0.27 


1.62 0.05 
2.64 0.15 
3.10 0.20 




bond angle FWHM 


bond angle FWHM 


0-Si-O 
Si-O-Si 


108° 12° 
141° 26° 


109.6° 10° 
142.0° 25° 



TABLE VL Statistics of hydrogen-bonding between the two surfaces and the first ML of water 
(within the first peak of the density distribution) . 



Surface 


OW-HW- • • OW 


OW-HW- • • OB OW HW- • • OH OH H- • • OW 


Total surface 


a-Si02 


2.54 


0.14 0.30 0.32 


0.76 


native 


2.51 


0.24 0.26 0.26 


0.76 


bulk water 


3.13 


- - - 


- 
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Figure Captions 

FIG. 1: Side and top views of tlie liydroxylated native oxide layer on Si(OOl) as obtained in a series 
of FPMD simulations.— li^ O atoms (red) either bridge two Si atoms or form part of a terminating 
hydroxyl group. The valence states of the 13 oxidised Si atoms range from Si^ to Si^^. 

FIG. 2: Examples of the Si-Si two-body interaction for Si'^^-Si'^^ species (dashed line) and Si^"*"- 
Si^+ species (dotted line), compared with the original SW interaction for uncharged Si (solid line). 
The equilibrium bond lengths are similar and both potentials tend to the Coulomb interaction for 
large separations. 

FIG. 3: The Si-0 two-body interaction. The equilibrium bond length increases with decreasing Si 
partial charge. 

FIG. 4: DFT (dashed lines) and classical (solid lines) binding energy curves of a single water 
molecule above the hydroxylated native oxide surface. In each simulation, the water molecule 
is displaced vertically, and the relative total energy is plotted as a function of its separation, d, 
from the surface. In the first configuration (left, black curves), the water molecule, with a H 
atom pointing downwards is bound to a bridging O atom. In the second configuration (right, 
blue curves), the water O atom is pointing downwards towards a Si atom and is repelled from the 
surface. 

FIG. 5: (left) Total radial distribution function in bulk amorphous silica. The first three peaks 
correspond to Si-0, 0-0 and Si-Si nearest neighbours, (right) 0-Si-O and Si-O-Si bond angle 
distributions, which are consistent with a network of corner-sharing Si02 tetrahedra. 
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FIG. 6: (left) Relative densities of fourfold coordinated (Si-4) and threefold coordinated (Si-3) Si 
atoms and bridging (0-2) and dangling (0-1) O atoms with increasing distance from the centre of 
the amorphous silica slab. Defects are concentrated in the first 5 A of the surface, (right) 0-Si-O 
bond angles near the surface (dashed line) show broadening out to larger angles than in the bulk 
(solid line), indicating the presence of Si-3 species. 



FIG. 7: Average density of water molecules perpendicular to the hydroxylated amorphous silica 
(left) and native oxide (right) surfaces. In both cases, the main peak is centred on the surface/water 
interface, with some water penetrating the surfaces and density oscillations continuing into the 
water layer. The water density in the centre of the simulation cells matches that of bulk water. 



FIG. 8: (top) Force-separation curves for two hydroxylated amorphous silica (left) and native oxide 
(right) surfaces. The bare surfaces (dashed lines) repel (positive force) at separations below A. 
Also plotted, from left to right, are curves corresponding to approximately 0.25, 0.50, 0.75, 1.00, 
1.50, 2.00 and 2.50 ML water coverage on each surface. The curves are translated upwards so that 
the force, at large separations, is zero, (bottom) Number of hydrogen bonds per water molecule 
as a function of surface separation for the amorphous silica (left) and native oxide (right) surfaces. 
Below about 1 ML (dotted lines), the number of hydrogen bonds does not reach the ideal surface 
density (upper dashed lines) before the surfaces repel. The maximum attractive force between 
surfaces at higher water coverages (solid lines) is limited by the hydrogen bond density between 
the surfaces, which tends to that of bulk water (lower dashed lines). 



FIG. 9: Surface energy of two hydroxylated amorphous silica (solid line) and native oxide (dashed 
line) surfaces as a function of water coverage. The bonding energy of silica is approximately 
constant, with a small peak at 1 ML, and decreases at low coverages. The native oxide surfaces are 
more strongly bound at low coverages, before tending to the silica and bulk water surface energies 
at high coverages. 
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FIG. 10: Average water density profiles perpendicular to two natively oxidised Si wafers at equi- 
librium separation, at increasing water coverage. The corresponding force-separation curves are 
numbered (top left). At 0.5 ML coverage (1), a large force is required to separate the single density 
peak into two. At 1.0 ML (2), the bonding energy is a maximum, because a stable water network is 
set up between the two wafers. At 2.5 ML (3), the interfacial water density and associated surface 
energy is close to that of bulk water. 
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